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Abstract 

A distributed model predictive control (DMPC) approach based on distributed 
optimization is applied to the power reference tracking problem of a hydro power 
valley (HPV) system. The applied optimization algorithm is based on accelerated 
gradient methods and achieves a convergence rate of O (p-) , where k is the itera- 
tion number. Major challenges in the control of the HPV include a nonlinear and 
large-scale model, nonsmoothness in the power-production functions, and a glob- 
ally coupled cost function that prevents distributed schemes to be applied directly. 
We propose a linearization and approximation approach that accommodates the 
proposed the DMPC framework and provides very similar performance compared 
to a centralized solution in simulations. The provided numerical studies also sug- 
gest that for the sparsely interconnected system at hand, the distributed algorithm 
we propose is faster than a centralized state-of-the-art solver such as CPLEX. 
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1. Introduction 

Hydro power plants generate electricity from potential energy and kinetic en- 
ergy of natural water, and often a number of power plants are placed along a long 
river or a water body system to generate the power at different stages. Currently, 
hydro power is one of the most important means of renewable power generation 
in the world In order to meet the world's electricity demand, hydro power 
production should continue to grow due to the increasing cost of fossil fuels. How- 
ever, hydro electricity, like any renewable energy, depends on the availability of 
a primary resource, in this case: water. The expected trend for future use of hy- 
dro power is to build small-scale plants that can generate electricity for a single 
community. Thus, an increasingly important objective of hydro power plants is 
to manage the available water resources efficiently, while following an optimal 
production profile with respect to changes in the electricity market, to maximize 
the long-term benefit of the plant. This water resource management must be com- 
patible with ship navigation and irrigation, and it must respect environmental and 
safety constraints on levels and flow rates in the lakes and the rivers. By sig- 
nificantly increasing the power efficiency of hydro power valley (HPV) systems, 
real-time control of water flows becomes an important ingredient in achieving this 
objective. 

An HPV may contain several rivers and lakes, spanning a wide geographical 
area and exhibiting complex dynamics. In order to tackle the plant-wide control of 
such a complex system, an HPV is often treated as a large-scale system consisting 
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of interacting subsystems. Large-scale system control has been an active research 
area that has resulted in a variety of control techniques, which can be classified in 
three main categories: decentralized control, distributed control, and centralized 
control. The application of these approaches can be found in a rich literature on 



control of water canals for irrigation and hydro systems [116 



14ll . We are inter- 



ested in applying model predictive control (MPC), a control method that has been 
successfully used in industry [|25|], thanks to its capability of handling hard con- 
straints and the simple way of incorporating an economical objective by means 
of an optimization problem. For the control problem of open water systems, 
centralized MPC has been studied in numerical examples using nonlinear MPC 
approaches in combination with model smoothing and/or model reduction tech- 



niques 1113 



systems [133 



! Ik and in real implementations with linear MPC of low-dimensional 



3411 . However, centralized MPC has a drawback when controlling 
large-scale systems due to limitations in communications and the computational 
burden. These issues fostered the studies of decentralized MPC and distributed 
MPC for large-scale water systems. Early decentralized MPC methods for irriga- 
tion canals used the decomposition-coordination approach to obtain decentralized 
versions of LQ control [6]. Several decentralized MPC simulations applied to 



irrigation canals and rivers were presented in UJl BO 
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261]. Distributed MPC 



approaches based on coordination and cooperation for water delivery canals were 



presented in [[7 



iflQQi. 



The typical control objective in these studies is to reg- 



ulate water levels and to deliver the required amount of water to the right place 
at some time in the future, i.e., the cost function does not have any special term 
except the quadratic penalties on the states and the inputs. On the other hand, in 
hydro power control, there are output penalty terms in the cost function that rep- 
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resent the objective of manipulating power production. Recent literature taking 
into account this cost function includes centralized nonlinear MPC with a paral- 
lel version of the multiple- shooting method for the optimal control problem using 
continuous nonlinear dynamics [29], and a software framework that formulates 
a discrete-time linear MPC controller with the possibility to integrate a nonlinear 
pred iction model and to use commercial solvers to solve the optimization problem 



i24|l . The hydro power control problem considered in the current paper is similar 



to the setup in ||2S 



2411 ■ However, it distinguishes itself by using a distributed 



control structure that aims to avoid global communications and that divides the 
computational tasks into local sub-tasks that are handled by subsystems, making 
the approach more suitable for scaling up to even more complicated hydro power 
plants. 

The distributed MPC design approach proposed in this paper is enabled by a 
distributed optimization algorithm that has recently been developed by the authors 
in [|8|] . This optimization algorithm is designed for a class of strongly convex prob- 
lems with mixed 1-norm and 2-norm terms in the cost function, which perfectly 
suits the power reference tracking objective in the HPV control benchmark. The 
underlying optimization algorithm in [8], although being implemented in a dis- 
tributed way, is proved to achieve the global optimum with an O(^) convergence 
rate, where k is the iteration number. This is a significant improvement compared 



to the distributed MPC methods presented in IS 141 121 USD, which achieve an O(^) 
convergence rate. There are three main challenges in applying distributed MPC 
using the algorithm from [8] to the HPV benchmark problem. The first one is that 
the nonlinear continuous-time model yields a relatively large linear model after 
spatial and temporal discretizations. We present a decentralized model order re- 
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duction method that significantly reduces the model complexity while maintaining 
prominent dynamics. The second challenge is that the power production functions 
are nonsmooth, which prevents gradient-based methods to be applied directly. A 
method to overcome this difficulty and to enable optimal control using the al- 
gorithm from [8] is also presented. The third challenge is that the whole system 
should follow a centralized power reference which, if the algorithm from [jsl] is ap- 
plied directly, requires centralized communication. We propose a dynamic power 
division approach that allows to track this centralized power reference with only 
distributed communications. By means of numerical examples, we will demon- 
strate the fast convergence property of the distributed algorithm which, when im- 
plemented on a single core, can outperform a state-of-the-art centralized solver 
(CPLEX) when solving the same optimization problem. 

The remaining parts of the paper are organized as follows. In Section |2l we 
describe the HPV system and the power reference tracking problem that were for- 
mulated in the HPV benchmark problem [|28ll . Section |3] provides a summary of 
the distributed optimization framework that the authors have developed in [80 . In 
SectionUl we present our approach for modeling and model reduction of the HPV 
system, followed by a reformulation of the MPC optimization problem, and de- 
veloping a distributed estimator so that the closed loop distributed MPC scheme 
can be implemented using neighbor-to-neighbor communications only. The sim- 
ulation results are presented in Section |5l which also features a comparison with 
centralized MPC and decentralized MPC. Through the various aspects of the com- 
parison including performance, computational efficiency, and communication re- 
quirements, the advantages of the distributed MPC algorithm will be highlighted. 
Section |6] concludes the paper and outlines future work. 
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2. Problem description 

In this section, we provide a summary of the hydro power valley benchmark 



[12 811 and we present the linearized model that serves as the starting point of our 



controller design. 

2.1. Hydro power valley system 

We consider a hydro power plant composed of several interconnected sub- 
systems, as illustrated in Figure [T] The plant can be divided into 8 subsystems, of 
which subsystem Si is composed of the lakes Li, L2, the duct Ui connecting them, 

-Ri, R2, respectively. Sub- 
system 5*2 is composed of the lake L3 and the ducts C2, T2 that connect L3 to the 
reaches -R4, -R5, respectively. There are 6 other subsystems, each of which consists 
of a reach and a dam at the end of the reach. These six reaches Ri, R2, R3, R^, R5, 
and Rq are connected in series, separated by the dams Di, D2, -D3, -D4, and D5. 
The large lake that follows the dam Dq is assumed to have a fixed water level, 
which will absorb all the discharge. The outside water flows enter the system at 
the upstream end of reach Ri and at the middle of reach R^. 

There are structures placed in the ducts and at the dams to control the flows. 
These are the turbines placed in the ducts Ti, T2 and at each dam for power pro- 
duction. In the ducts Ci, C2 there are composite structures that can either function 
as pumps (for transporting water to the lakes) or as turbines (when water is drained 
from the lakes). 

The whole system has 10 manipulated variables, which are composed of six 
dam flows (qoi, qD2, Qds, Qoi, Qob, Qdo), two turbine flows {qxi, Qtt) and two 



'a reach is a river segment between two dams. 
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pump/turbine flows (qci, ^02)- Further, the system has 9 measured variables, the 
water levels in the three lakes {hn, hL2, hiz) and the water levels at the end of 
each reach {hnx, /ih2, /^ra, hm, Hr^, hm)- 



One of the control problems specified in [28'] is the power reference tracking 
problem. We introduce state variables x, which consist of water levels in the lakes 
and reaches and water flows within the reaches, and control variables q, which are 
the manipulated water flows. The problem is to track a power production profile, 
p^'^^{t), on a daily basis using the following cost function: 




Figure 1 : Overview of the HD-MPC hydro power valley system 0281] 



2.2. Power reference tracking problem 
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8 J' 

+ E / (^^(^) - - €)dt (1) 

subject to the nonlinear dynamics and linear constraints on outputs and inputs as 
specified in [12 811 . The weights Qj, -Rj, i = 1, . . . , 8, 7, and the testing period T are 
parameters of the benchmark. 

The quadratic term in the cost function represents the penalties on the state 
deviation from the steady state x^^ and the energy used for manipulating the in- 
puts away from the steady state flows q^^. The 1-norm term represents the power 
reference tracking mismatch, in which the function p^'^^ is the power reference and 
the function represents the locally produced/consumed power by a subsystem 
i G {1, . . . , 8}. For 2 = 1,2 the produced/consumed power is (cf. [28]) 

where qc^ and qx^ are the flows through ducts Ci and Tj, Axc, and Axt^ are the 
relative differences in water levels before and after ducts Ci and respectively, 
is the power coefficient of the turbine Tj, and 

(3) 

is a discontinuous power coefficient that depends on whether the duct Ci acts as a 
turbine (qc^ (t) > 0) or as a pump (gc-. (t) < 0). For i = 3, . . . , 8 we have 

Pi{x{t),q{t)) = kD,_2qD,.2{t)^XD,_2{t) (4) 

which is the power produced by the turbine located at dam -Dj 2- The pro- 
duced/consumed power functions given in ^ and dH) are nonlinear, and even 
nonsmooth for subsystems 1 and 2 due to the differences of k^c. and kp^ in (|3]), 
thus complicating a direct application of a standard MPC scheme. 
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Still, the complexity of the system and control objective suggests an optimiza- 
tion based control strategy, such as MPC. Further, the distributed nature of the 
system makes it possible to consider distributed MPC techniques. However, the 
stated optimization problem ([T]) is a nonlinear continuous-time dynamic optimiza- 
tion problem, which in general is very hard to solve. In the next sections we will 
discuss the modeling of the hydro power valley that leads to a linearized model. 

2.3. Nonlinear hydro power valley model 

The model of the reaches is based on the one-dimensional Saint Venant partial 
differential equation, representing the mass and momentum balance (see [28] for 
details): 

{ dqjt, z) ^ ds{t, z) ^ ^ 
ld'fq{t,z)\\ 1 d (q\t,z)\ dh{t,z) 

with z the spatial variable, t the time variable, q the river flow (or discharge), 
s the cross-section surface of the river, h the water level w.r.t. the river bed, /f 
the friction slope, Io{z) the river bed slope, and g the gravitational acceleration 
constant. 

The partial differential equation (|5]) is converted into a system of ordinary 
differential equations by using spatial discretization. To achieve this, each reach 
is divided into 20 cells, yielding 20 additional states, which are the water levels at 
the beginning of the cells. For details of the spatial discretization and the equations 
for the resulting nonlinear dynamical system the reader is referred to [28, Section 
2.1.1]. The resulting nonlinear dynamical system has in total 249 states, 10 inputs, 
and 9 outputs. 
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2.4. Model linearization and discretization 

As mentioned in Section [23] a set of nonlinear ordinary differential equations 



that describe the hydro power valley dynamics is presented in [1281 Section 2.1.1]. 



A linear continuous-time model which is linearized around the steady state oper- 



ating point (x*^*^, q^^) is also provided in the HPV benchmark package [|28ll . Dis- 
cretizing this model using zero-order-hold gives a discrete-time linear system with 
249 states and 10 inputs. The coupling of the subsystems is through the inputs 
only. This implies that discretization using zero-order-hold of the continuous-time 
system keeps the structure of the original system description. Thus, the resulting 
discrete time system has a block-diagonal dynamics matrix, a block-diagonal out- 
put matrix, and a sparse input matrix, and each subsystem i = 1, . . . , 8 can be 
expressed in the following form: 



{k + l)= A,,x1{k) + %??(^) (6) 

in which the variables x'^,q^, and stand for the deviation from the steady- state 
values, and the subscripts i,j stand for the subsystem indices. As mentioned 
the subsystems are coupled through the inputs only and at least for some j G 
{1, . . . , 8} we have Bij = for every i = 1, . . . ,8. 

The use of a discrete-time linearized model enables controller design with 
some specific approaches, which include our proposed distributed optimization 
technique presented in [|81]. Before describing our main contributions, we now 
provide a summary of this distributed optimization framework in the next section. 
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3. Distributed optimization framework for MPC 

In this section, we describe the distributed optimization algorithm developed 
in [|8D which is based on an accelerated gradient method. The first accelerated 
gradient method was developed in [I21I1 and further elaborated and extended in 



23i 



320. The main idea of the algorithm presented in 180 is to exploit the 
problem structure of the dual problem such that accelerated gradient computations 
can be distributed to subsystems. Hence, the distributed algorithm effectively 
solves the centralized optimization problem. Dual decomposition has been used 
in the past to tackle the complexity of large-scale optimization problems arising 
in water supply networks pi]. In our work however, in addition to simplifying the 
local computations, we apply this decomposition philosophy in order to distribute 
the decision-making process. 

The algorithm in [|81] is developed to handle optimization problems of the form 

min ^x^Hx + 5(^x + 7||xa||i (7) 

X,Xa I 

s.t. Ax = b 
Cx< d 

Xa = PX - p 

where x G M" and Xa G are vectors of decision variables, and x is partitioned 
according to: 



r T 

[xi , . . . , X 



MJ 5 



(8) 



and Xj G M"'. Further, the matrix H G M"^" is positive definite and block- 
diagonal, the matrices A G M'^^", C G M''^", and P G M™^" have sparse struc- 
tures, and g G W\ p G W\ b G M^, d G W . We introduce the partitions g = 
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[gl, glY, p = [pf , . . . , vlY, b = K , . . . , hlV, d 



H 



Hi 



Cii 



Hm 

CiM 



C 



A/1 



c 



MM 



'M\ ' 

An 

Ami 
Pii 

Pa/i 



FjT jTIT 
[til ) • • • ; ^AfJ ' 

Aim 

• Ama/" 

PlAf 

• PAfAf 



where the partitions are introduced in accordance with (|8} and Qi E R"' , Pi G M*"' , 

The assumption on sparsity of A, C and P is that Aij = 0, Cjj = 0, and Pij = 
for some i, j and we assume that the constraint matrices are built such that An ^ 
0, Cjj 7^ 0, and ^ for alH = 1, . . . , M. Based on the coupling, we define 
for each subsystem a neighborhood set, denoted by Mi, as follows: 

Mi = {j G {1, . . . , M}| A,j 7^ or Aji 7^ or ^ or Cji ^ or (9) 

P,,- T^OorP,, ^0}. 

Note that there are two type of equality constriants in (|7), the first one involves 
only X and the matrix A has a sparsity pattern, i.e., there is no global coupling 
introduced in that equality constraint; the last one involves both x and Xa, more- 
over introduces a global coupling due to the fact that Xa is penalized in the 1-norm 
term of the cost function, thus it is not straightforward to deal with this constraint 
as we could treat the first constraint. Throughout the paper, the dual variables 
corresponding to these constraints are treated differently, and a distributed ap- 
proximation of the 1-norm term is introduced to treat the second type of equality 
constraint. 
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We introduce dual variables A G W, n E W,!/ E for the equality con- 
straints, inequality constraints, and equality constraints originating from the 1- 
norm cost in <^} respectively. We also introduce the dual variable partitions 
A = [Af , . . . , XIj]'^, fi = [/if, . . . , nL]'^, and u = [uf, where A.; E M«% 
Hi E W\ and z/j E M™'. Based on [|8D, the dual problem of dV]) can be cast as the 
minimization of the negative dual function 

/(A, /i, v) = i(A^A + + P^vf-R-'iA^X + CV + P'^^)+ (10) 

+ b'^A + + V^i^ 

and the dual variables are constrained to satisfy 

XeW, /xGR;, z/g[-7,7]'" (11) 



where IR+ denotes the non-negative real orthant. The negative dual function (iTOl) 
has a Lipschitz continuous gradient with constant (cf. [8]) 

L = ||[A^ P^]^H-i[A^ P^]||2 (12) 

and can hence be minimized using accelerated gradient methods. The distributed 
accelerated gradient method as presented in dsl] is summarized below in a slightly 
different form that is adapted to our HPV application problem at hand. 

Algorithm 1. Distributed accelerated gradient algorithm 

Initialize X^ = X^^, fi^ = fi^^, = u^^ and x^^ with the last values from the 
previous sampling step. For the first sampling step, these variables are initialized 
by zeros. 

In every node, i, the following computations are performed: 
For A; = 0,1,2,... 
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1. Compute 



= - 



xf = xf + J ^(xf - xf ^) 



2. 5en<i xf to each j G A/i, receive y.^ from each j G A/i 



3. Compute 



Af + |^(A?-A-) + i(5:A,.,-b,) 

/^r' = max jo, /.f + |^(/if - /.f-i) + ^ ( E - d.) | 



z/f^+i = mill < 7, max 



+ i(gP,.,-p,)]} 



4. Af +\ /if +\ z/f +' to eac/z j G receive A^+\ /i^+\ z/^^+' /rom each 
J e Afi. 

The Lipschitz constant L of V/ is used in the algorithm. For MPC purposes 
we only need to compute L once in a centralized way and use it through all MPC 
problem instances. 

Besides the suitability for distributed implementation, another merit of Algo- 
rithm [His its fast convergence rate. The main convergence results of Algorithm [T] 
are given in [8], stating that both the dual function value and the primal variables 
converge towards their respective optima with the rate of O (^) where k is the 
iteration index. This convergence rate is much better than the convergence rate of 
classical gradient-based optimization algorithms, which is O (^). 
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4. Control of HPV using distributed MPC 

We have so far described the linear discrete-time model of the HPV in Sec- 
tion |2] and the fast distributed optimization method, Algorithm [H that serves as a 
basis for designing a distributed model predictive controller to be applied to the 
HPV. However, there are three major challenges for this application. First, the lin- 
ear discrete-time model cannot be directly used in an MPC context due to the ex- 
istence of a number of unobservable and uncontrollable modes. These unobserv- 
able/uncontroUable modes are a result of the spatial discretization in each reach 
which creates states that cannot be observed/controlled separately. In addition, the 
linear discrete-time model has a large number of states, causing a large computa- 
tional burden. Second, the power functions associated with the ducts Ci and C2 
are nonsmooth (cf. © and ©). The nonsmoothness is caused by the fact that the 
flow through Ci and C2 is bidirectional and the powers consumed/produced do 
not have equivalent coefficients. The third major challenge is the global coupling 
in the cost function due to the fact that we have to track a central power refer- 
ence function that specifies the desired sum of locally generated power outputs. 
This global coupling prevents a distributed implementation of Algorithm [T] since 
the sparsity in the constraints is lost. These issues are addressed in the following 
sections. 

4.1. Modification of the linear model 

In this section we show how to create a model of the HPV that is suitable 
for the DMPC framework presented in [8]. First we present a model reduction 
technique that keeps the system structure, then the nonsmooth power function is 
treated. 
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4.1.1. Decentralized model order reduction 

The block-diagonal structure of discrete-time dynamical system Q makes it 
possible to perform model reduction on each subsystem individually. Several 



model reduction methods have been proposed for interconnected systems 



35, 



27D 



In this work, we use a straightforward balanced truncation method [IllL 1 1711 to re- 
duce the order of each local model (|6]). 

Let us introduce Bi = [Bn . . . Bis] and q"^ = [{qf)'^ . . . (gf )^]^ to get the 
following discrete-time linear model of each subsystem: 

xf{k + l) = Auxf{k) + Biq''{k) (13) 
yf{k) = axfik)- 

Applying the balanced truncation technique yields transformation matrices de- 
noted by T/ and T/'^"^ for each subsystem, where j^rj^r.mv _ j denoting the 
new state variables, xj = T[xf, and the control variable q^ = q^, we represent the 
reduced order model as: 

xl{k + l)=Alxl{k) + Blq^{k) (14) 
ym = Clxlik) (15) 

where AJ, = T^'AnT-''''^, B] = T^Bi and Q = C.t;'''''^. It should be noted that 
the block- sparsity structure of Bl is the same as in the non-reduced input matrix 
Bi, since the model reduction is performed for each local model separately. More- 
over, all the modes of the reduced model are both observable and controllable. 

The model reduction gives a 32-state reduced model that approximately rep- 
resents the dynamics of the full linear model with 249 states. 
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4.1.2. Treatment of nonlinear and nonsmooth power function 

One of the difficulties in applying a linear MPC approach to the hydro power 
valley is the nonsmoothness of the power function associated with the ducts Ci 
and C2, which is included in the expression for power generation (|2]) in subsys- 
tem 1 and subsystem 2, respectively. In order to handle this nonsmoothness, we 
use a double-flow technique, which means introducing two nonnegative positive 
variables to express the flow in Cj, z = 1, 2 at a sampling step k: 

• qcip{k): virtual flow such that Ci functions as a pump 

• QCiT (^)' virtual flow such that Ci functions as a turbine 

The introduction of virtual flows requires the input-matrices, Bj, to be augmented 
with two extra columns identical to the ones multiplying qci,i = 1,2 with the 
opposite sign to capture that pump action is also introduced with a positive flow. 
The resulting reduced order model has 12 inputs instead of the original 10. Using 
the introduced flows and qc„, the power function (|2]) for subsystems 1 and 2 
can be rewritten as 

Pi{x{k),q{k)) = {kTc^qc„{k) - kp^^qc,p{k)^ ^Xc^k) + kr^qT^k) Axr^k) 

(16) 

with the additional constraints that g^xl^) ^ ^^QCtpik) > and qc„{k)qcip{k) = 
0. The last constraint expresses the fact that water flows in only one direction 
at a time, i.e., that either the pump or the turbine is active. The resulting non- 
linear expression (fT6] ) can in turn be linearized around the steady- state solution 
(x^^,q^^). Since q^_ = for i = 1,2 we get the following linear local power 
production/consumption approximation for subsystems i = 1,2: 



17 



Pi{x{k),q{k)) = Ax\ 



+ 



+ kT.q'r] (A^T,(fc) - AxfJ + %Ax- {grXk) - g??) + 

+ kr^q^Ax'^^ 

This reformulation results in a linear expression with a nonlinear constraint at 
each time step k, qc„{k)qcip{k) = 0, that approximates the original nonsmooth 
nonlinear power production/consumption expression (|2j. We show our approach 
to handle the nonlinear constraint in Section 14.21 

For subsystems i = 3, . . . , 8 we have smooth power production expressions 
(01) that can be directly linearized without introducing virtual flows: 

Pi{x{k),q{k)) = kDSDAx% + kDS% (Axd.(A;) - Axl^J + 

+ kDA^% {qDAk)-qD,) 

4.2. HPV optimization problem formulation 

In this section we formulate an optimization problem of the form (|7]) that can 
be used for power reference tracking in the HPV benchmark using MPC. We have 
obtained a linear discrete-time dynamical system (fT4])- (fT5l) for the HPV with state 
variables x^ and control variables q^ . The constraints are upper and lower bounds 
on the outputs and inputs and their values can be found in UM . Using the trans- 
formations matrices T/ and T/'"™, these constraints can readily be recast as linear 
constraints for the reduced order model variables x\q^ . The power reference 
problem formulation ([T) specifies a quadratic cost on states and control variables 
and a 1-norm penalty on deviations from the provided power reference, p^'^^. For 
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control horizon, A^, this optimization problem can be written as 

N-i r 8 



t=o I i=l 



s.t. fll), ([15]) k = 0,...,N -1 i = l,...,8 

Clxlik) k = 0,...,N-l i = 1, . . . , 8 

qi{k)eQ^ k = 0,...,N-l t = l,...,8 

x,ik) = (A;) - Zti M^'ik), q'{k)) k = 0,...,N-l 
qc,Ak)qcAk) = k = 0,...,N-l ^ = 1,...,2 
where and Qi are sets representing the local output and input constraints, the 
additional variable Xa captures the power reference tracking mismatch, and the 
notation x represents the stack of variables x- (A;) and ql{k) for all i and k, while Xa 
is the stacked variable of Xa(/i;) for all k. Note that we can write x = [xf , . . . , x|^]^ 
where each Xj, i = 1, . . . , 8 includes all the variables that belong to subsystem i. 

4.2.1. Power reference division 

Since the original cost function contains a non-separable 1-norm term, the 
power reference constraints in the optimization problem dTTl ) are coupled between 
all subsystems. This implies that Algorithm [T] requires some global communica- 
tion even though the only information needed to be sent to the global coordinator 
is pi{x^{k), q^{k)) for A; = 0, . . . , — 1 from each subsystem z = 1, . . . , 8. 

In order to obtain a suitable dual problem, we first need to reformulate the 
cost function in a separable form. For the sake of brevity, we focus on one sam- 
pling step and drop the time index k. Thus for now our simplified objective is to 
decompose the following problem: 



mm 

{Xi}i=l,...,8 



p'-^^-^P.x 

i=l 
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(18) 



with X = [x^, . . . , Xg and Pi the matrix coefficient such that the power function 
produced or consumed by each subsystem pi(a;'^(/c), q^{k)) is linearized as Pix{k). 

In this section we present two different ways that avoid global communica- 
tion when solving this problem. In the first approach, we divide and distribute 
the global power reference to the subsystems in a static fashion. In the second 
approach, we show how the subsystems can trade local power references between 
neighbors to achieve a satisfactory centralized reference tracking. 

Static local power references. The idea here is straightforward. We divide the 
global power reference into local ones, i.e., p^''^ is divided into local parts pj^^ 
i = 1, . . . , 8. We have chosen to compute pf^ such that it satisfies 

f^'] = 3^-^"^^-^^^^ , for^ = 0,...,iV-l (19) 

with Pi{x^^, q^^) the power produced by subsystem i in the steady-state condition. 

This means that the fraction of the total power reference given to subsystem i 
is constant. The optimization problem is changed accordingly, i.e., the following 
cost function can be used instead of (fTSi) : 

8 



min 2, 

1=1 



(20) 



with X = [xj^, . . . , Xg ] . This allows for a distributed implementation since the 
matrix Pj introduces only local couplings, i.e., subsystem i needs only neighbor- 
ing and local water levels and local water flows to compute the corresponding 
power output. The disadvantage of the static power reference division is that the 
global power reference tracking is not very accurate, as will be shown in the sim- 
ulations section. 
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Dynamic local power references. The static power division essentially means that 
each subsystem always tracks a fraction of power reference that is equal to the pro- 
portion it produces in the steady-state condition. When the total power reference 
deviates significantly from the steady-state power, this idea may not work well 
since the proportional change of the local power reference can lead to sub-optimal 
performance. Inspired by an idea in [15], we now introduce the dynamic power 
division, in which the subsystems have more flexibility in choosing the appropri- 
ate local power reference to be tracked. The main idea is that each subsystem will 
exchange power references with its direct neighbors. 

Let us define for each pair with j G M a node that is in charge of 
determining the power exchange variable between subsystems i and j, denoted 
by 5ij if node i is in charge and by 5ji if node j is in charge 1^. Then for each 
subsystem we form the set^: 



with 6i the vector containing all j G Aj, and the nominal power reference 
for subsystem i. In words, the local power reference for each subsystem i devi- 
ates from the nominal value by adding the exchange amounts of the links that i 

^Note that here we discuss the power division for each sampHng step, i.e., there are 5ij (fc) or 

5ji{k) withfc = 0,...,A^-l. 

simple way is to let the subsystem with smaller index lead the exchange, i.e., Ai — G 



= {j I J ^ is in charge of 5ij}. 



(21) 



Now we replace (fT8l) by the following cost function: 




8 



(22) 



N^,j > i}. 



21 



manages and subtracting the exchange amounts of the links that affect i but are 
decided upon by its neighbors. Note that problem (|22] ) has a sparse structure that 
complies with the existing sparse structure of the HPV system, i.e., this method 
does not expand the neighborhood set of each subsystem. 

The advantage of this dynamic power division is that it makes use of the exist- 
ing network topology to form a sparse cost function, and the total power reference 
is preserved even if the local power references can deviate from the nominal val- 
ues, i.e., we always have: 



Now that we have a separable cost function by using either a static or a dy- 
namic power division technique, we can cast the approximate optimization prob- 
lem in the form (|7]) that has a separable dual problem, and apply Algorithm \T\ at 
every sampling step. However, due to the requirement of positive definiteness of 
the quadratic term in the objective function, the introduced power exchange vari- 
ables Sij must be penalized with a positive definite quadratic term. This implies 
that power reference exchange has an associated cost. 

Communication structures. In the preceding sections we have presented three dif- 
ferent ways to handle the power reference term. The first is the one with cen- 
tralized power reference term which we hereby denote by GLOBAL-REF. The 
second is the one with static local power references which we denote by LOC- 
REF-STAT. The third is the dynamic local power reference which from here on is 
denoted by LOC-REF-DYN. In Table[l]we provide an overview of the neighbor- 
hood sets Mi for the different power reference tracking schemes. 




(23) 
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Table 1: Neighborhoods of subsystems (A/i) 



Subsystem 


GLOBAL-REF 


LOC-REF-DYN 


LOC-REF-STAT 


1 


{!,■ 


.,8} 


{1,3,4} 


{1,3,4} 


2 


{!,• 


.,8} 


{2,6,7} 


{2,6,7} 


3 


{!,■ 


.,8} 


{3,1,4} 


{3,1,4} 


4 


{!,• 


.,8} 


{4,1,3,5} 


{4,1,3,5} 


5 


{!,- 


.,8} 


{5,4,6} 


{5,4,6} 


6 


{!,■ 


.,8} 


{6,2,7,5} 


{6,2,7,5} 


7 


{!,- 


.,8} 


{7,2,6,8} 


{7,2,6,8} 


8 


{!,■ 


.,8} 


{8,7} 


{8,7} 



We can see that all subsystems have the same neighborhood sets for the dy- 
namic local reference tracking and the static local reference tracking. 

4.2.2. Relaxation of nonlinear constraint 

The second issue that hinders the optimization problem (fTTT ) from being solved 
using Algorithm [Hare the nonlinear constraints qci.^{k)qc^p{k) = with i = 1,2. 
In this section we present a way to relax these constraints. 

Assuming in the cost function we have the penalty Rd [qc,t1c,p]'^ on the pump 
and turbine action in ducts Ci, i = 1,2, with 



R, 



Rci 



T 

R, 







■Cii 



(24) 



We also have the constraints that qcipik) > 0, qcn{k) > and qc,T{k)qCip{k) = 
0. We relax this by removing the nonlinear constraint and adding a cross-penalty 
a^Rc^pRciT for some a G (0, 1) in the cost function, i.e., we set 



R. 



Ra- 



(^\/Rc,pRciT 



a 



^jRc^pR. 



R. 



(25) 
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This relaxation is implementable using the proposed algorithm since the nonlinear 
constraint is removed and replaced by a cross-penalty. The cross-penalty gives an 
additional cost if both qc„ and qc^p are non-zero. The closer a is to 1, the larger 
the penalty. For a > 1 it is easily verified that we lose strong convexity on the 
quadratic cost function, i.e., Rd loses positive definiteness and such choices for 
a are therefore prohibited. 

The relaxation is not equivalent to the original nonlinear constraint and thus 
cannot guarantee that the nonlinear constraint is respected using this relaxation. 
However, it turns out that the optimal solution using the cross-penalty in the cost 
(l25l) in most simulated cases coincides with the optimal solution when the nonlin- 
ear constraint qc„ {k)qci:p (k) = and the original diagonal cost (|24] ) are enforced. 
In some cases however, the optimal solution using the relaxation does not respect 
the nonlinear constraint. To address this, a two-phase optimization strategy is 
developed and presented next. 

4.2.3. Two-phase optimization 

We propose a two-phase optimization strategy as an ad-hoc branch and bound 
optimization routine that uses two consecutive optimizations. In the first opti- 
mization the relaxed optimization problem is solved. If the nonlinear constraints 
are respected, i.e., we get a solution that satisfies gc'.^(A;)gc'.p(A;) = 0, the global 
optimal solution for the non-relaxed problem is found. If some of the nonlin- 
ear constraints do not hold, the optimization routine is restarted with setting the 
smaller flow between qc„ (k) and qc^p (k) to zero, for i = 1,2, k = 0, . . . , N — 1. 
The resulting algorithm is summarized below. 



Algorithm 2. Distributed branch and bound algorithm 
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1 . Solve the relaxed problem using Algorithm\T\ 

2. Mqc„{k)qcAk)^^, z = 1, 2, t = 0, . . . , iV - 1 

liqc^k) > qcM 

Add constraint: qc^p{k) = 

Else 

Add constraint: qc^^{k) = 

End 

End 

3. Solve relaxed problem using AlgorithmUjwith the additional flow constraints 

This ad-hoc branch and bound technique does not theoretically guarantee that 
the optimal flow directions are chosen. However, we can guarantee that the non- 
linear constraints are always satisfied. Further, for the distributed MPC formu- 
lation we will see in the simulations section that the global optimal solution for 
the non-relaxed problem is found at every time step using this branch and bound 
algorithm. 

4.3. Distributed estimation 

From Section |2] we know that not all states can be measured, which implies 
that an observer needs to be used to feed an initial condition to the optimizer. 
The reduced-order linear model (fT4T i-(|T5t has local dynamics and outputs only, 
which implies that an observer can be designed in decentralized fashion. We 
introduce the local estimate and the local observer-gain K^, and the following 
local observer dynamics 

xlik + 1) = Al,xm + Blq^-{k) + KM{k) - Clx^,{k)) 
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Because of the sparse structure of Bl this observer can be implemented in a dis- 
tributed fashion where only the inflows to subsystem i need to be communicated. 
The estimation error = x] — x\ has local error dynamics 

^{k + l) = {Al-K,Cl)x\{k) 

Thus, the observer can be designed in a decentralized fashion and be implemented 
in a distributed fashion. 



5. Simulation results 

We perform distributed MPC simulations of the hydro power valley using 
3 different ways of handling the power reference: GLOBAL-REF, LOC-REF- 
DYN, and LOC-REF-STAT, using the proposed Algorithm |2] We also solve the 
problem ([17]) using a state-of-the-art MIQP- solver, namely CPLEX. In CPLEX 
the nonlinear constraints given in (flTT ) can be addressed by introducing binary 
variables. More specifically, for each duct Ci,i = 1,2, we define two virtual 
flows, qa^p and qc,T, and require that both values are nonnegative. Each virtual 
flow has a maximum capacity, hence the constraints for these flows are: 

< < q^:^ ^^^^ 

< gc„ < q^:^ 

We introduce binary variables bi E {0, 1} and impose the following constraints: 

qc-,<qd:^{l-b.) 

The constraints (|26l ) and (|27T ) ensure that either g^-p = 0, gc-^ > (if = 1) 

or qc„ = 0, qc,p > (if bi = 0). 
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This formulation results in an MIQP for which there are efficient Branch-and- 
Bound algorithms implemented in CPLEX. To make the 1-norm term in (fTTT ) fit 
the MIQP-formulation used in CPLEX we introduce auxiliary variables v and use 
the following equivalent reformulation 

min||Pa; — -x^ minify 

X x,v 

S.t. — V < Px — p < V 

We also compare the proposed distributed MPC method to a decentralized MPC 
approach in which each subsystem solves its own local MPC problem without 
any communication, in order to show the advantage of DMPC w.r.t. decentralized 
MPC. 

5.1. Simulation details 



We use the original nonlinear continuous model presented in [|28ll as simula- 
tion model. The ode-solver odelSs in MATLAB is used to perform the simula- 
tions. A MATLAB function that computes the derivatives needed by odelSs is 
provided in the benchmark package [28J. The control system consists of the dis- 
tributed observer from Section |43l which feeds Algorithmic] with estimates of the 
current state. 

Besides the mismatch between the model used for control and the model used 
for simulation we have also added bounded process noise to capture mismatch 
between the simulation model and the real plant. The magnitude of the worst 
case process noise was chosen to be 1% of the steady-state level x^^. We also use 
bounded additive measurement noise where the measured water levels are within 
±3 cm from the actual water levels. 
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(c) DMPC and LOC-REF-DYN 




(d) DMPC and GLOBAL-REF 



Figure 2: Comparison of power reference tracking performance using DMPC and decentralized 
MPC approaches. Solid Unes: produced power, dashed lines: reference power, dotted lines: steady 
state power. 



We use a sampling time of 30 minutes in all simulations and the control hori- 
zon is = 10, i.e., 5 hours. The simulations are performed over a 24 hour period 
since the power reference trajectories are periodic with this interval. 

All simulations and optimizations were implemented on a PC running MAT- 
LAB on Linux with an Intel(R) Core(TM) 17 CPU running at 3 GHz and with 4 
GB RAM. 
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5.2. Control performance comparison 



The power reference tracking results are plotted in Figures |2(d)f|2(a)] where 
the full power reference and the sum of the local power productions are plot- 
ted. The scheme GLOBAL-REF achieves very good tracking performance, while 
the scheme LOC-REF-STAT shows a significant deterioration in tracking perfor- 
mance. However, the introduction of the possibility to exchange power references 
in LOC-REF-DYN between subsystems restores the very good tracking perfor- 
mance while keeping the computations distributed. The tracking performance of 
the decentralized MPC approach is very poor, due to the lack of communications. 
Hence, it is recommended not to use a decentralized MPC approach, unless com- 
munication is prohibited due to the lack of communication facilities or due to the 
policy of different authorities. 



In Appendix A and Appendix B there are figures that show the input and 
output evolutions and the corresponding constraints with the scheme LOC-REF- 
DYN. We can observe that all constraints are satisfied despite disturbances, model 
mismatch, and the use of an observer. For the schemes GLOBAL-REF and LOC- 
REF-STAT all the constraints on the inputs and outputs are also satisfied. 

During the simulations, it is observed that all schemes achieve stable closed- 
loop behaviours, which can be explained that the HPV system is already marginally 
stable and does not have critical dynamics, and the prediction horizon is long 
enough so that the MPC controllers do not introduce instability to the closed loop. 
Note that neither the centralized MPC nor the distributed or decentralized MPC 
approaches used in this simulations employ a method that provides guaranteed 
stability to the closed-loop system, since this property is beyond the scope of this 
paper. Based on the techniques for distributing the computation and improving 
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the efficiency of the algorithm that are proposed in this paper, one can further 
incorporate other MPC schemes that guarantee the closed-loop stability, which 
could be important for other types of applications where there are large mismatch 
between the nonlinear and the linearized models. 

5.3. Computational efficiency/accuracy 

In Table [21 we provide a comparison of the execution times of the central- 
ized MPC problems (flTl) . We compare the distributed Algorithm [2] to the solver 
CPLEX when solving (fTTl ). i.e., with power-division GLOBAL-REF in Algo- 
rithm 121 To solve this problem using CPLEX, an MIQP formulation is used. In 
every iteration of Algorithm |2l the relaxed problem is solved twice. We also com- 
pare the above execution times to the case when we solve the first relaxed prob- 
lem in Algorithm|2l which is a QP, using CPLEX. At each sampling step, the same 
problem is solved, and the execution time t is measured. Although in this example 
the solvers easily solve the problem within the time frame of the sampling time, 
we can see that the computation time for our MATLAB -implemented algorithm is 
always lower than the C-implemented CPLEX for both the MIQP and QP cases. 

Table 2: Comparison of computation time between Algorithm|2]and CPLEX for 48 instance of the 
same problem 





Algorithm |2l 


CPLEX for MIQP 


CPLEX for QP 


min t (s) 


0.023 


0.087 


0.049 


max t (s) 


0.086 


0.121 


0.089 


average t (s) 


0.054 


0.098 


0.063 


std dev t (s) 


0.017 


0.009 


0.009 



As previously discussed. Algorithm |2l cannot guarantee that the global op- 
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timum for (flTT ) is found. However, in the DMPC simulations presented in this 
section the global optimum of (fTT]) is found at every sampling step using Algo- 
rithm |2l 

5.4. Communication requirements 

The sizes of the optimization problems using power reference division GLOBAL- 
REF, LOC-REF-DYN or LOC-REF-STAT are almost equal. Comparing GLOBAL- 
REF to LOC-REF-STAT we get some additional constraints due to the power 
reference division and comparing LOC-REF-DYN to LOC-REF-STAT we get 
some additional decision variables 5ij to enable distributed power reference re- 
assignment. 

In Table [3] the number of iterations njter needed to obtain the solution is pre- 
sented. The average and max values of riiter and the standard deviation are com- 
puted using 48 simulation steps, i.e., 24 hours. 



Table 3: Number of iterations to solve the MPC optimization in one step 





Alg.[I]with 
GLOBAL-REF 


Alg.[I]with 
LOC-REF-DYN 


Alg.[I]with 
LOC-REF-STAT 


average niter 


311.3 


579.1 


942.5 


max riiter 


498 


1054 


2751 


std dev niter 


93.8 


210.9 


440.8 



We can notice that different DMPC schemes converge with different average 
numbers of iterations. The reason is that for LOC-REF-STAT it is more difficult 
to satisfy the different 1-norm terms with equality, i.e., to follow the local power 
references. This implies that the corresponding dual variable u becomes large 
(close or equal to 7) and it takes more iterations to achieve convergence. As 
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a result, the scheme LOC-REF-STAT with a simpler communication structure 
might require more communication resources than e.g., GLOBAL-REF, which 
has a more complicated communication structure but needs fewer iterations. 

In order to estimate the total time required for communications within each 
sampling time, we now assume the worst case happens in every iteration of Algo- 
rithm m in which Algorithm [T] has to be executed two times. In Algorithm \T\ 
also assume the worst case that every primal and dual variable has to be ex- 
changed between distributed controllers, with prediction horizon = 10 there 
are 10 x (44 + 65) = 1090 variables to be transmitted once per iteration. Let each 
variable be a 32-bit floating-point, then the total time it would take for transmitting 
exchanged variables in 1000 iterations is: 

2 X 1090 X 32 X 1000 = 69, 760, OOO(bits) (28) 

or roughly 70 Mbits. With a decent wireless network that can connect each two 
nodes with a rated transfer as 7 Mbps, the total time for communications is less 
than 10 seconds for one thousand iterations. Note that in practice, there should be 
more communication delays due to the initialization of transmissions. Since the 
communication time is considerably shorter than the sampling time of 30 minutes, 
the iterative methods taking about one thousand iterations sampling time can still 
be implemented in real time. 

The scheme LOC-REF-DYN performs very well in terms of communication, 
computation, as well as performance aspects and is therefore the chosen candidate 
for distributed implementation for the given case study. 
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6. Conclusions and future work 

The proposed distributed MPC approach has been applied to the power ref- 
erence tracking problem of the HD-MPC hydro power valley benchmark. Two 
distributed schemes have been compared to centralized and decentralized MPC 
methods. We have provided relaxations and approximations for the original non- 
linear nonsmooth problem formulation as well as proposed a way to follow a cen- 
tralized power reference in a distributed fashion. Furthermore, we have presented 
a practical branch-and-bound algorithm that solves all optimization problems en- 
countered in the simulations and achieves as good performance as the centralized 
MPC that is known to have global optimum. The simulation results show that 
the introduced approximations and relaxations capture the behavior of the system 
well and that very good control performance is achieved. Finally, a comparison 
to state-of-the-art optimization software (CPLEX) shows that the proposed algo- 
rithm has significantly better execution times in general. 

As the next step before implementation in real plants, the proposed distributed 
MPC approach should be tested against different hydraulic scenarios and other 
HPV setups. To cope with varying water flows entering the system, these should 
be estimated and compensated for. Furthermore, a weather model could be in- 
cluded that estimates the future inflows to the system. 
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Figure A. 3: Input constraint satisfaction using Algorithm |2] and power division LOC-REF-DYN. 
Dash-dotted lines: upper bounds, dashed lines: lower bounds. 
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Figure B.4: Output constraint satisfaction using Algorithm|2]and power division LOC-REF-DYN. 
Dash-dotted lines: upper bounds, dashed lines: lower bounds. 
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